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Overview 

In this paper, we describe a computational fluid dynamics (CFD) technique based on unstructured 
triangular/tetrahedral meshes. A finite-volume scheme is used in conjunction with a multi-stage 
Runge-Kutta algorithm. Convergence enhancements in the form of dual time-stepping and time- 
derivative preconditioning are used to overcome the limitations of conventional multi-stage schemes. 
The method is applied to propulsion-related flows and shown to perform satisfactorily. 


Technical Discussion 

Governing Equations and Discretization 

The equations governing the unsteady flow of a multi-species fluid mixture may be written in 
integral Cartesian form for an arbitrary control volume, V y with surface, S, as follows: 

w,IJJ WdV + JJ FdS+ JI GdS = JJJ Hdv <» 

where 

W = [p, pv x , pv y , pv z , pe, pYif 

F = [/9®, pv x v + pi, pvyv + pj, pv z v + pit, pev + pv, puVj] 

G = [ 0 , r xi , Tyi, T zi , TijVj + q , pv t Yi] T 

H = [0, 0, 0, 0, ui] T 

Here p y v, e, and p are the density, velocity, total energy, and pressure of the fluid mixture, 
respectively. Y\ and u)/ are the mass fraction and rate of production of the 1-th chemical species, r, 
q and vi are the viscous stress tensor, and the heat flux and diffusion velocity vectors, respectively. 

The domain is divided into triangular/ tetrahedral volumes and the governing equations dis- 
cretized over these volumes using a cell-centered finite-volume approach. 5 Odd-even decoupling due 
to the use of central-difference operators is damped by adding fourth difference artificial dissipation; 2 
a first-order dissipative term is selectively added near discontinuities. 

Preconditioning for Convergence Enhancement 

The performance of the time-marching scheme at low Mach numbers is enhanced by using time- 
derivative preconditioning 1, 4 * 7 in the context of unstructured meshes. Equation (1) is modified as 
follows: 

4 l JiI Wdv+ JI Fds+ JJ° ds =JJI Hdv 
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where W are a new set of primary dependent variables: W - [p,v T ,v y ,v z ,T,Yi] T , and T is the 
preconditioning matrix: 
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m 
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0 0 0 0 0 

p 0 0 0 0 

0 P 0 0 0 

0 0 P 0 0 

pv x pVy pv z pCp p{ht - h N ) 

0 0 0 0 P 


Here U T is a reference velocity, U T < Min(\v\,U c ), where U c is a characteristic speed of the flow 
such as the local speed of sound or a maximum velocity within the domain. 

In the preconditioned system, time-step definition is based on the modified eigenvalues: A±, |w|, M,.~ 

where 

A ± = l[H(l + K/3t/ r 2 )±c'] 

C ' = V /b| 2 ( 1 -«^r) 2 +4f/ r 2 
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dp It pCp 


The parameter k is included to provide time-step control at low Reynolds numbers. 1 When vis- 
cous effects become important, the preconditioning scheme alters the acoustic speed such that the 
CFL number is of the same order of magnitude as the von-Neumann number; thus optimal wave 
propagation speeds as well as optimum von-Neumann numbers result. 7 


Dual Time-Stepping for Unsteady Flows 

To provide for efficient, time accurate solution of the governing equations, we employ dual time- 
stepping, 3 adapted for use with an explicit multi-stage scheme. Here we introduce a preconditioned 
pseudo- time derivative term into (1) as follows: 


*IJ J J* dV + liJ J J wdV+ II FdS+ J J GdS = J 1 J Hiv (2) 

Note that as r -»• oo, the first term on the LHS of (2) vanishes and (1) is recovered. 

The time-dependent term in (2) is discretized in an implicit fashion by means of a second order 
accurate, three-point backwards difference in time and the pseudo-time derivative is driven to zero 
by means of the following multi-stage algorithm: 


r + ( w" 

2 At dW 


w (0) = 

w w 

-W (0) ) = 

—ajAr jltW 

W {k+l) = 



= -a i Ar|RW ( ‘" 1) + ^(3W (,_1) -4W ( " ) + iy (n 1) )| 


( 3 ) 


Here t = 1,2 
pseudo-time 


, ..., m is the stage counter for the m-stage scheme, and k and n represent any given 
and physical-time level, respectively. Throughout the iterations on k, W and 


136 



W^ n ^ are held constant. As r — » oo, — ► JV^ n+1 \ Note that the matrix on the LHS 

of (3) is inverted in a point- wise fashion and its inverse is readily derived analytically and need not 
be computed numerically. Note also that the physical time step, A t, is limited only by the level 
of desired temporal accuracy. And the pseudo time-step, A r, is determined by the multi-stage 
scheme, which here employs local time stepping and preconditioning for convergence enhancement. 

Results 

Rocket Engine Flowfield 

To demonstrate the viability of using unstructured solution-adaptive meshes to compute transonic, 
internal, viscous flows typically found in space propulsion applications, the finite-volume, multi- 
stage method presented here is used to compute a transonic, converging/diverging nozzle flowfield. 
Figure 1 shows the unstructured triangular mesh used in these computations. This mesh consists of 
6,877 cells. It has been adapted to gradients of velocity in order to resolve the boundary layer along 
the nozzle walls, and to gradients of pressure in order to better capture the oblique shock in the 
divergent section. Pressure contours within the nozzle, plotted on a logarithmic scale in Figure 2, 
depict the oblique shock that develops as the expanding flow is turned inward by the nozzle walls. 

Flow Past Circular Cylinder 

The time accurate, dual time-stepping scheme described above is applied to solve the unsteady, two- 
dimensional flow over a circular cylinder in crossflow at a Reynolds number, R e = Ud/v = 75. At 
this Reynolds number the flow is laminar and exhibits periodic unsteady behavior as vortices shed 
from the cylinder to form a Karman vortex street in the wake. The computational domain for this 
problem is chosen to extend 5 diameters upstream and 20 diameters downstream of the cylinder, 
with symmetry boundaries placed 5 diameters above and below. Details of the unstructured mesh 
in the vicinity of the cylinder are shown in Figure 3(a). The specified fluid is air at standard 
conditions with U = 11.5 m/s (M = 3xl0“ 4 ) and d = 1 cm. The predicted shedding frequency, /, 
is 1.67s” 1 , resulting in a Strouhal number, St — fd/U = 0.146. This compares reasonably with 
the measured value 6 of 0.147. Contours of stream function in the vicinity of the cylinder at several 
times during the shedding cycle are show in Figure 3(b). The physical time step, At is chosen to 
be 0.025s and time periodic behavior is achieved in roughly 5s beginning from an impulsive start 
from rest at t = 0. Twenty-five inner iterations are performed at each time level to achieve three 
orders of magnitude decrease in the x-momentum residual using the preconditioned, multi-stage 
scheme with a CFL of 2.7. Thus 600 iterations are required to resolve one time period of 0.6s. 
This represents a 1000 time speed-up over global time stepping which, for this grid and these flow 
conditions, would have restricted the time step in each cell to roughly lxl0 -6 s. 

Conclusions 

A finite-volume, multi-stage algorithm based on an unstructured grid topology is presented. The 
use of solution adaption is demonstrated in the calculation of a typical space propulsion application. 
Dual time-stepping and time-derivative preconditioning are shown to provide efficient solution of 
unsteady, low speed flow on an unstructured mesh. The benefits of extending this method to low 
Mach number flows with heat release, and to incompressible flows with variable density are evident, 
and our efforts continue in that direction. 
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Figure 1: Unstructured, triangular mesh used for transonic, converging/diverging nozzle calcula- 
tions. 
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Figure 2: Pressure contours in convergent /divergent nozzle (plotted on logarithmic scale) computed 
using unstructured, triangular mesh. 
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